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EXECUTIVE  SUMMARY 


Military  operations  in  cold  climates  require  that  soldiers  be  maintained  relatively  warm  for 
safety  and  optimum  performance.  Protection  afforded  by  passive  cold  weather  clothing 
ensembles  may  be  limited  in  duration.  This  problem  is  accentuated  under  inactive  conditions 
and,  particularly,  in  relation  to  the  extremities  and  the  exposed  parts  of  the  face.  These  body 
elements,  even  when  protected  by  thermal  insulation,  e.g.,  gloves,  may  cool  down  to 
temperature  levels  which  may  endanger  the  soldiers’  safety  and  ability  to  operate.  Knowledge 
of  endurance  times  associated  with  these  exposures,  is,  therefore,  of  prime  importance  to  the 
military  commander.  In  this  study,  the  estimation  of  endurance  times,  as  it  pertains  to  the 
digits,  is  performed  by  a  one-dimensional,  cylindrical  model  depicting  these  body  elements. 

An  analytical  solution  of  the  axial  temperature  variations  is  obtained  in  terms  of  an  infinite 
series.  Thermal  and  physical  properties  specified  for  the  model  pertain  to  the  modeled  body 
parts  and  to  typical  insulation  layers  applied.  Blood  perfusion  effects  are  lumped  into  a 
volumetric  heat  generation  term.  Cold-induced  vasodilatation  (CIVD)  effects  were  not 
included  in  the  present  analysis.  Endurance  times,  as  determined  by  a  drop  of  cylinder  (digit) 
tip  temperature  to  5°  C  were  evaluated.  Parameters  included  in  this  evaluation  were:  (a) 
environmental  temperatures,  (b)  thermal  insulation  applied  on  the  cylinder,  (c)  length  of  the 
cylinder,  and,  (d)  diameter  of  the  cylinder.  It  was  found  that  the  lower  the  temperature,  the 
longer  the  finger,  and  the  smaller  the  diameter  thereof  -  the  shorter  is  the  endurance  time  for 
the  same  thermal  insulation.  Results  of  the  model  were  compared  to  measured  data  for  a 
subject  not  exhibiting  CIVD  response  to  cold  stress.  Conformity  of  calculated  and  measured 
data  was  very  good  with  a  maximum  deviation  of  less  than  10%  at  only  one  particular  point 
in  time.  This  model  facilitates  the  conservative  estimation  of  lower  bounds  to  thermally 
insulated  fingers  and  toes  exposed  to  cold  stress. 
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INTRODUCTION 


Humans  are  often  required  to  occupy  and  operate  in  cold  climates.  The  limited  natural 
ability  to  cope  with  extreme  environmental  conditions  must  be  supplemented  by  auxiliary 
means  if  humans  are  to  survive  and  function  effectively  under  harsh  conditions.  One  of  the 
ccm-'cn  means  employed  to  cope  with  cold  environments,  is  the  use  of  clothing  to  increase 
r  .rr  ?„  insulation.  !n  many  cases  the  protection  arfordsd  by  clothing  is  sufficient  but  there 
may  be  situations  wherein  certain  body  parts,  most  notably  the  extremities  and  the  face,  are 
still  not  adequately  protected.  This  may  occur  particularly  under  prolonged  sedentary 
conditions  in  very  cold  environments  [1].  This  report  analyzes  the  temperature  changes  of 
the  extremities,  particularly  fingers  (digits)  and  toes  during  cold  stress. 

The  reasons  for  the  susceptibility  of  these  body  elements  to  cold  temperatures  are  both 
physiological  and  physical  in  nature.  The  physiological  reason  relates  to  the  main  heat 
source  in  an  extremity,  e.g.,  finger,  which  is  its  blood  supply.  Under  neutral  environmental 
conditions,  the  flow  of  blood  and  the  associated  internal  generation  of  metabolic  heat  are 
usually  sufficient  to  maintain  finger  temperatures  within  a  safe  and  comfortable  range.  Upon 
exposure  to  cold  temperatures,  however,  blood  flow  to  the  fingers  is  constricted,  primarily  as 
a  means  for  conserving  heat.  This  causes  a  significant  diminution  of  the  supply  of  heat  to 
the  finger.  Values  in  the  literature  indicate  a  drop  in  blood  perfusion  rate  from  about  15-40  to 
only  about  1  cc/(min  100  cc)  [2  cited  in  3].  An  even  higher  range  of  variability,  of  about  100 
fold,  was  also  suggested  [1], 

One  of  the  physical  reasons  for  the  susceptibility  of  fingers  to  cold  environments  is  due  to 
their  shapes,  which  resemble  slender  cylinders.  These  cylinders  act  like  heat  transfer  fins. 
The  result  is  enhanced  heat  loss  to  the  environment,  with  an  ensuing  accelerated  drop  in 
finger  temperature.  This  physical  property  of  an  extremity  is  not  specific  to  cold  exposures. 
Its  effects,  however,  are  much  more  pronounced  under  such  conditions,  mainly  due  to  the 
increased  temperature  difference. 

An  additional  physical  limitation  relates  to  the  thickness  of  thermal  insulation  which  may 
be  applied  on  cylindrical  surfaces.  Since  the  surface  of  a  cylinder  is  curved,  each  layer  of 
applied  insulation  would  be  of  a  larger  diameter  and  hence  of  an  increased  surface  area. 
Thus,  there  is  a  "critical"  thickness  beyond  which  heat  loss  to  the  environment  would  be 
increased,  rather  than  decreased.  This  limitation  does  not  apply  to  insulation  applied  to  flat 
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surfaces  due  to  the  absence  of  curvature.  In  the  context  of  this  study,  gloves  may  be 
considered  to  provide  individual  layers  of  insulation  to  each  finger.  Mittens  are  different  in 
that  they  usually  cover  all  four  fingers  at  once,  except  the  thumb.  They  can,  therefore,  be 
made  thicker,  or  more  insulative,  compared  to  the  size  limitations  imposed  on  the  finger 
covers  in  gloves. 

Van  Dilla  et  al.  [3]  have  estimated  the  insulation  required  to  maintain  finger  temperatures 
at  comfortable  levels  of  15.5°C  (60°F)  or  28.9°C  (84°F)  under  different  internal  heat 
generation  rates  and  environmental  temperatures  down  to  -40°C  (-40°F).  They  performed  a 
steady  state  analysis  and  simulated  the  application  of  fabrics  to  either  cylinders,  spheres  or 
plane  surfaces,  depicting,  finger  shaft,  fingertip  and  hand,  respectively.  Assuming  that  "one 
quarter  of  an  inch  of  fabric  thickness  may  be  taken  as  the  maximum  thickness  consistent 
with  reasonable  dexterity",  the  authors  concluded  that  the  finger  could  be  protected  by  a 
maximum  insulation  value  of  about  0.9  do  [3].  This  finding  demonstrates  the  limitations 
imposed  by  physical  laws,  on  the  one  hand,  and  the  physiological  response  of 
vasoconstriction,  on  the  other.  Van  Dilla  et  al.  further  extended  their  analysis  to  mittens 
concluding  that  "mittens  are  decidedly  superior  to  gloves  in  both  insulation  value  and 
efficiency  of  fabric  utilization".  This  conclusion  holds,  however,  at  the  expense  of  rendering 
the  mittens  "highly  cumbersome,  unwieldy  and  impractical". 

Using  Hatch’s  equation  of  body  cooling,  Van  Dilla  et  al.  estimated  the  tolerance  times  for 
the  average  finger  skin  temperature  to  drop  to  8.9°C  (48°F)  from  an  initial  temperature  of 
33°C  (91.4°F),  while  covered  by  a  certain  mitten.  They  found  that  at  -40°C  ambient 
temperature,  the  tolerance  times,  for  different  heat  generation  rates,  were  80-85  min 
compared  to  experimentally  determined  75  min.  At  -23.3°C  (-10°F)  and  -6.7°C  (20°F)  the 
calculated  tolerance  times  were  156-184  min  and  132  min,  respectively.  These  calculated 
times  compared  favorably  with  those  determined  experimentally,  180  min  and  137  min., 
respectively,  for  these  two  conditions. 

Molnar  analyzed  the  rate  of  digital  cooling  to  the  freezing  point  [4].  He  assumed 
Newton’s  law  of  cooling  to  apply  and  a  uniform  temperature  throughout  the  digit.  Fingers 
were  exposed,  without  insulation  in  a  small  "wind  tunnel”,  to  a  -15°C  environment  with  wind 
at  10  m/sec.  Finger  temperatures  were  observed  to  fall  to  -6.2°C,  or  even  lower,  before 
freezing  occurred  followed  by  a  quick  rewarming  to  about  -2.4°C  due  to  the  "evolution  of 
(latent)  heat".  In  certain  experiments,  conducted  under  similar  conditions,  no  freezing  was 
observed  to  have  occurred,  indicating  the  onset  of  cold  induced  vasodilatation  (CIVD). 
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lAcTiar  concluded  his  study  by  stating  that  "there  was  no  correlation  between  the  elapsed 
ume  to  me  onset  of  freezing  and  the  relative  cooling  rale",  he  also  suggested  the  possibility 
ihat  CIVD  started  when  a  particular  skin  temperature  was  reached  but  refrained  from  giving 
any  further  details. 

Wilson  and  Goldman  [5]  and  Wilson  et  al.  [6]  studied  the  freezing  temperatures  of  finger 
skin.  Their  results  indicate  that  in  the  majority  of  cases,  50  out  of  a  total  of  72,  skin  freezing 
,io  occur  under  the  conditions  of  the  study.  The  mean  supercooled  skin  temperature  at 
wf,;„h  frostnip  appeared  was  -9.4°C  followed  by  an  average  temperature  rise  of  5.3°C  due  to 
the  absorption  of  the  heat  of  fusion  at  ice  crystallization  [6j.  Wilson  and  Goldman  suggest 
that  below  a  wind  chill  of  1400  and  for  air  temperatures  in  the  range  of  -10°  to  -15°C, 
freezing  of  fingers  rarely  occurs  [5]. 

Goldman  discussed  the  problem  of  protecting  the  inactive  Arctic  soldier  for  8  hours  in  a  - 
40°C  environment  with  a  1.34  m/sec  (3  mph)  wind  [1].  As  a  consequence  of 
vasoconstriction,  the  hands  and  feet  become  the  limiting  factors  in  tolerance  to  cold.  He 
showed  that,  under  these  conditions  no  passive  protection  was  possible  and  that  auxiliary 
heating  was  required.  Results  show  that  electrical  resistance  heating  of  gloves  and  socks 
may  be  the  most  practical  means  of  maintaining  the  extremities  at  any  desired  temperature 
level  between  10°C  (50°F)  and  40.2°C  (105°F).  The  minimal  power  requirements,  for 
maintaining  the  extremities  above  4.4°C  (40°F),  were  3  W  and  7  W  for  the  hand  and  foot, 
respectively.  Goldman  found  that  the  mean  times  to  cool  from  maintained  5th  finger 
temperatures  of  26.7°C  (80°F)  or  15.6°C  (60°F)  to  4.4°C  (40°F),  after  power  was  shut  off  to 
the  heated  glove,  were  6.7  and  4.1  min,  respectively.  When  a  windproof  leather  glove  shell 
was  provided  over  the  heated  glove,  these  times  became  slightly  longer  at  9.6  and  5.9  min, 
respectively. 

An  alternative  to  auxiliary  heating  of  the  hands  and  feet  in  cold  environment  is  metabolic 
heat  generated  by  physical  activity.  This,  however,  may  not  always  be  totally  effective, 
depending  on  the  environmental  conditions  and  level  of  activity.  Santee  and  co-workers  [7], 
studied  the  effects  of  cold  weather  handwear  used  by  petroleum  handlers  who  were  required 
to  operate  manual  pumps  intermittently.  Their  experiments  were  conducted  in  a  -34.4°C  (- 
30°F)  and  a  2.2  m/s  (5  mph)  wind  speed  environment.  Initial  results  indicated  that 
endurance  times,  determined  when  a  finger  temperature  reached  5°C,  barely  exceeded  10 
min.  When  both  pre-test  conditions  were  adjusted  and  test  chamber  temperature  was  raised 
lo  -28.9°C  (-20°F),  endurance  times,  including  pumping  times,  increased  considerably  to 
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between  about  27  and  33  minutes  The  autho  's  showed  that  middle  finger  nail  bed 
temperature  vaned  opreciably  in  subjects  demonstrating  the  fast  initial  drop  followed  by 
oscillations  due  to  the  alternating  pump-rest  cycles. 

The  purpose  of  the  present  report  is  to  analyze  the  factors  affecting  the  thermal  response 
of  a  fmger-like,  cylindrical  model.  A  one-dimensional  analytical  model  of  the  axial 
temperat|,re  /ariations  of  this  model  is  formulated.  The  model  is  solved  by  an  orthogonal 
transformation  and  is  presented  in  terms  of  an  infinite  series.  The  effects  of  finger  insulation, 
length  and  diameter  and  those  of  the  environmental  temperature,  were  studied.  Analytical 
predictions  of  the  model  were  compared  to  experimental  data  and  a  very  good  conformity 
was  obtained. 


ANALYSIS 


Consider  the  right  angle  circular  cylinder,  shown  in  Fig.  1,  to  depict  any  body  element 
which  may  be  approximated  by  this  geometry,  e.g.,  fingers,  toes,  arms,  etc.  This  cylinder  is 
exchanging  energy  with  the  environment  through  its  exposed  surfaces.  Inside  the  cylinder 
heat  flows  from  one  location  to  another  only  by  conduction.  In  addition,  there  is  heat 
exchange  at  the  base  of  the  cylinder  due  to  its  attachment  to  the  proximal  body  element, 
e.g.,  hand  in  the  case  of  a  finger. 

For  simplicity,  the  following  assumptions  are  made: 

(a)  The  thermophysical  properties  of  the  cylinder  are  uniform  and  isotropic. 

(b)  Metabolic  heat  is  generated  at  a  uniform  rate  in  the  cylinder,  which  may  vary  in  time. 

(c)  Due  to  the  aspect  ratio  of  the  cylinder,  radial  heat  conduction  is  ignored  as  compared 
to  axial  conduction. 

(d)  Blood  perfusion  effects  on  heat  exchange  are  not  included  directly  in  the  model,  but 
are  lumped  into  the  term  describing  metabolic  heat  generation. 

Conservation  of  energy  requires  that  the  following  balance  be  maintained  in  each 
infinitesimal  element  of  the  cylinder: 

[Cl LANGE  1N\  [  ENERGY  |  f  ENERGY  )  f  ENERGY  \  [ ENERGY  LOST \ 

{  STORED  \  =  {CONDUCTED}  -  {CONDUCTED)  *  {GENERATED}  TO  THE  > 

(  ENERGY  j  (  INTO  J  (  OUT  OF  j  [  INSIDE  j  [ENVIRONMENT) 
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In  mathematical  notations  this  heat  balance  may  be  written  as: 


P 


K*l 

dz2 


t0) 


(1) 


where:  p  is  density,  c  is  specific  heat,  T  is  temperature,  t  is  time,  k  is  thermal  conductivity  ,  z 
is  axial  coordinate,  q  is  metabolic  heat  generation  rate  per  unit  volume,  p  and  A  are  cylinder 
circumference  and  cross-section  area,  respectively,  h  is  heat  transfer  coefficient  between  the 
cylindrical  surface  and  the  environment  and  T0  is  environmental  temperature. 

Equation  (1)  is  solved  subject  to  the  following  boundary  and  initial  conditions: 

At  the  base  of  the  cylinder  a  variable  temperature  is  assumed 

7=^(0  §z  =  0,/i0  (2) 


At  the  tip  of  the  cylinder,  convective  heat  exchange  with  the  environment  occurs: 

=  h,(T-  T0)  @  z  =  «.  f*0  (3) 


where  h,  is  the  heat  transfer  coefficient  between  the  tip  of  the  cylinder  and  the  environment 
which  may  be  different  than  that  assumed  for  the  circumferential  area  of  the  cylinder,  h. 

The  initial  condition  for  the  axial  temperature  in  the  cylinder  is  given  by  an  arbitrarily 
specified  distribution 


T  =  T,{z)  @  t  =  0,  0  <;  z  s  C 


(4) 


The  solution  to  the  above  set  of  equations  is  obtained  by  employing  an  orthogonal 
integral  transformation  [8]  and  [9].  Steps  performed  in  the  derivation  of  the  solution  are 
discussed  in  Appendix  A.  The  result  obtained  is 
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(5) 


T{z,  l ) 


om  n  ^ 

-  a  f  G(r\)Qxp{akzy])dr\}exp(-ak2t)s\n~^— 

„.i  J0  ( 

t0  -  m  -  t0\ 


wnate, 


(6) 


2  (P„2  +  Bi2) 
^(P  n^Bi2  +  Bi) 


(7) 


G(0  - 


^<LF0(t) 

P/7 


1 

a 


d*b}  +  <7(0  «*, 

dt  +  *Pn 


(cos  p„  -  1) 


(8) 


Fo(0  =  r,w  -  r0 


(9) 


'■«  ’  ifs/'f)  - 1  (,0) 

and, 

P„cotp„  =  -Bi  (11) 

The  general  solution,  Eq.  (5),  contains  a  number  of  operations  to  be  carried  out.  These 
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operations  are  either  integrals  in  Eqs.  (5)  and  (6),  or  the  solution  of  a  transcendental 
eigenvalue  equation  (10).  Performance  of  all  indicated  operations  depends  on  the  specific 
value  of  B,  in  Eq.  (10),  the  specific  initial  condition  T,  (z),  Eq.  (6),  and  the  values  specified  for 
both  the  time  dependent  boundary  condition  at  the  base  of  the  cylinder,  T,  (t),  and  the 
metabolic  rate,  q  (t),  in  Eq.  (8).  These  operations  would  now  be  carried  out  for  a  set  of 
selected  values  and  functions.  Obviously,  different  values  and  functions  may  be  selected 
and  substituted  into  the  appropriate  expressions. 

Solution  of  the  eigenvalue  equation  (10),  may  be  obtained  by  a  number  of  well- 
established  techniques.  Two  parameters  will  have  to  be  specified:  B,  and  the  acceptable  cut¬ 
off  accuracy  to  be  employed  by  the  solution  algorithm.  In  the  present  case,  a  combination  of 
solution  techniques  was  chosen,  each  one  is  applied  for  a  different  range  of  the  parameter 
B,.  For  low  and  negative  values  of  B|<1 .37,  a  Newton-Raphson  procedure  is  employed.  For 
B,>1 .37,  a  parameter-adjusted  axis-crossing  routine  is  applied.  In  all  cases,  care  is  taken  to 
insure  that  all  successive  eigenvalues  are  calculated  and  the  cut-off  accuracy  chosen  is 
usually  of  the  order  of  10  4,  particularly  for  the  first  ten  eigenvalues.  The  calculation 
procedure,  employing  double-precision,  was  programmed  in  Turbo-Pascal  and  is  included  in 
the  listings  presented  in  Appendix  B.  A  sample  output  of  results  is  given  in  Table  1. 

Next,  the  integration  indicated  in  Eq.  (6),  is  carried  out  for  a  specific  example: 


T 


rjsin—dir  - 


-  T„) 

r  n 


BJ 


—  f  T,(z)  sin  +  T0 cosp„  -  T:(o) 
Pn  o 


(12) 


It  is  assumed  that  the  initial  temperature  distribution  in  the  cylinder  is  linear: 
r,(z)  =  (T3  -  Tjj  ♦  Tt  (13) 
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which  yields,  upon  substitution  and  integration: 


T 


T^|-5nr^sinp"  *  r*  -  r'  -  <r3  -  T-o) cos p„) 

Pn  P  n 


(14) 


Equation  (14)  holds  also  for  the  specific  case  of  an  initially  uniform  temperature  by 
assigning  T2=T3.  Other  forms,  including  piecewise  initial  temperature  distributions,  may  be 
accommodated. 

Finally,  time-dependent,  exponentially-decaying  functions  are  assumed  for  the  cylinder 
base  temperature  T,(t)  and  the  metabolic  heat  generation  rate  q(t),  respectively: 

T,(t)  =  T,  +  (7,  -  T,)exp(-f/a)  (15) 

<7(0  =  qt  +  {q,  -  qf)exp{-tle)  (16) 


where  subscripts  f  and  i  indicate  final  and  initial  values,  respectively,  and  5  and  e  are  the  time 
constants  for  the  base  temperature  and  metabolic  heat  generation  rate,  respectively. 

In  summary,  and  for  the  specific  functions  [Eqs.  (13),  (15)  and  (16)]  chosen  here,  the 
temperature  of  the  cylinder  is  given  by 


T(z,()  -  r0  .  f  £  <(*„2X  ♦  C„  ♦  G„  ♦  H„]exp(-aa„20 

77-1 

sinP»7 

-  C„  -  G„exp(  f/0)  *  H„8XP( -(/.)) - r' 

-  Ijfgjj  -  1l(r'  *  (r<  -  r,)exp(-(f«)  -  r„) 
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TABLE  1:  Sample  array  of  the  first  40  eigenvalues  for  Bi=0.6 


INDEX 

EIGENVALUE,  pn 

1 

1 .868875893336337 

2 

4.830676037164286 

3 

7.926291752979685 

4 

11.047499501381431 

5 

14.177642504858346 

6 

17.311913171169223 

7 

20.448423440435814 

8 

23.586283165340159 

9 

26.725021611627494 

10 

29.864355859361979 

11 

33.004119522373991 

12 

36.144200808375331 

13 

39.284523617986736 

14 

42.425034356955287 

15 

45.565694191642798 

16 

48.706474295923925 

17 

51 .847352823254987 

18 

54.988312915192573 

19 

58.129341354865495 

20 

61.270427634230991 

21 

64.411563294015238 

22 

67.552741447665682 

23 

70.693956432130392 

24 

73.835203547724092 

25 

76.976478861650733 

26 

80.117779057723041 

27 

83.259101320085364 

28 

86.400443242290720 

29 

89.541802755508613 

30 

92.683178071327191 

31 

95.824567635801671 

32 

98.965970092251780 

33 

102.107384250924127 

34 

105.248809064086231 

35 

108.390243605450252 

36 

111.531687053073284 

37 

114.673138675067463 

38 

117.814597817595924 

39 

120.956063894738790 

40 

124.089269053266150 

RESIDUAL  OF:  pn  cot(pn)-Bi 

-1 .82260505531 735566E-0005 
7.6221 5668461 966775E-0005 
1 .053051 34443524871 E-0005 
2.790936961 62852848E-0006 
1 .02861 7451 92203398E-0006 
4.623096531 83909462E-00C7 
2.371 49491 793284348E-0007 
1 .33650781 39831 4024E-0007 
-8.75384552623320253E-0005 
-6.973361 7078901 0436E-0005 
-5.66851 208 1 0870751 0E-0005 
-4.680783871 06063377E-0005 
-3.91 231 265649041 849E-0005 
-3.300031 2912923751 3E-0005 
-2.801 81 5467081 591 1 0E-0005 
-2.388661 3384058261 9E-0005 
-2.040071 1764531 8699E-0005 
-1 .741 2265968761 51 60E-0005 
-1 .481 1 999561 3488296E-0005 
-1 .251 7899395786391 8E-0005 
-1 .04674449592920748E-0005 
-8.61 230079254881 892E-0006 
-6.91 462260674405565E-0006 
-5.34444082346335284E-0006 
-3.877771 4987731 2774E-0006 
-2.49523848554927038E-0006 
-1.181 04864801 584679E-0006 
7.7771 54857504901 82E-0008 
1.291 96504065475382E-0006 
2.4704792333291 4922E-0006 
3.62082767459327324E-0006 
4.7493431 270071 31 86E-0006 
5.861 4050931 2643243E-0006 
6.961 599557761 71 763E-0006 
8.053859845361 52299E-0006 
9.14156880406797332E-0006 
1 .022765857603791 07E-0005 
1.131 46708342631 462E-0005 
1 .24048257638330025E-0005 
1 .025944360891 801 08E+0000 
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where, 


*  |(cosp„  -  1)) 


(18) 


BJ  4(r,  -  T,){*6L  -  1) 
P  n  aX2n  6  -  1 


(IS) 


Wn 


[cos  Pn  -  1) 


lU(Qt  ~  Qr) 

aX2e  -  1 


(20) 


and  x,  X2,  Pn,  Bi  and  Bn  are  defined  by  Eqs.  (12),  (A.32),  (11),  (A.16)  and  (7),  respectively. 
Equation  (17)  was  programmed  in  Turbo  Pascal  and  a  listing  is  given  in  Appendix  B. 


RESULTS  AND  DISCUSSION 


The  present  model  may  be  studied  by  selecting  a  set  of  parameters  and  evaluating  its 
behavior.  The  parameters  should  be  selected  to  depict,  as  closely  as  possible,  realistic 
values  pertaining  to  the  modeled  entity.  In  the  present  case  these  parameters  should 
represent  biological  tissues  and  their  thermophysical  as  well  as  anthropometric  properties.  In 
addition,  there  are  external  conditions  to  be  specified,  e.g.,  environmental  temperature, 
insulation  applied,  model  initial  condition,  etc. 

The  number  of  possible  combinations  of  parameters  is  almost  infinite.  The  advantage 
afforded  by  an  analytical  model  is  in  its  ability  to  facilitate  the  study  of  specific  parameters  at 
a  relatively  low  cost  and  without  the  need  for  elaborate  and  costly  experiments.  However, 
the  ultimate  test  of  a  model  is  in  its  ability  to  closely  predict  the  behavior  of  the  modeled 
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entity.  This  may  be  tested  by  comparing  model  predictions  to  experimental  results. 

In  the  present  study,  three  cases  were  studied  and  are  presented  below.  All  cases  were 
assumed  to  relate  to  the  behavior  of  fingers,  i.e.,  digits,  covered  by  gloves  of  different 
insulating  values.  Parameters  employed  in  the  calculations  are  listed  in  Table  2. 

In  the  first  case,  a  constant  digit  base  temperature  and  metabolic  heat  generation  rate  are 
assumed.  Results  are  shown  in  Figs.  2  and  3,  for  an  environmental  temperature  of  -5°C. 
Digit  base  temperature  is  maintained  at  30°C  throughout  the  exposure  and  a  total  heat 
transfer  coefficient  of  7.12  W/m2  °C  representing  an  insulation  value  of  about  0.9  clo.  The 
initial  temperature  in  the  digit  is  assumed  to  be  linear  from  30°C  to  20°C,  as  shown  in  Fig.  2. 
The  other  curves  in  this  figure  show  the  temporal  evolution  of  the  axial  temperature 
distributions  along  the  digit  until  a  final  steady  state  is  reached.  As  is  to  be  expected,  digit 
temperatures  decrease  monotonically  until  finally  a  large  temperature  gradient,  of  about  30°C 
in  this  case,  has  been  established  between  the  base  and  the  tip  of  the  digit. 


TABLE  2:  Typical  parameters  used  in  the  calculations 
(parentheses  indicate  default  values) 


Thermal  conductivity,  W  m 1  C 1 

0.418 

Thermal  diffusivity,  m2  hr'1 

0.0004546 

Fleat  transfer  coefficient,  W  m 2  C 1 

5  -  12 

Digit  length,  m 

0.06  -  0.12  (0.08) 

Digit  diameter,  m 

Metabolic  heat  generation  rate,  W  m 3 

0.01  -  0.0175  (0.015) 

initial,  or  constant 

15000  -  29000  (15000) 

final 

15000  -  26500  (15000) 

Environmental  temperature,  C 

0  to  -20 

Initial  digit  temperature  (linear),  C 
Exponential  time  constants,  hr 

30  to  20 

metabolism 

0.4 

digit  base  temperature 

1.3 
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Figure  3  shows  temperature  variations  at  different  locations  along  the  digit.  The 
uppermost  curve  is  plotted  for  Z=0,  or  digit  base,  whereas  the  lowest  curve  indicates  digit  tip. 
The  total  time  for  which  both  Figs.  2  and  3  were  calculated  is  about  3.05  hrs.  At  this  time  a 
final  steady  state  has  been  established  throughout  the  entire  digit.  However,  inspection  of 
Fig.  3  reveals  that,  for  all  intended  purposes,  a  practical  steady  state  becomes  established 
much  sooner,  at  about  1  hour  and  30  minutes.  This  is  also  indicated  in  Fig.  2  by  the 
convergence  of  the  curves  at  the  final  temperature.  Thus,  the  curve  indicating  the  final 
temperature  distribution  in  the  digit  at  3.05  hrs  is  quite  close  to  that  calculated  for  less  than 
an  hour  (t/t*=0.3)  from  the  beginning  of  the  exposure. 

One  of  the  more  practical  results  that  may  be  deduced  from  Fig.  3  relates  to  the  time  it 
would  take  to  reach  a  certain  temperature  anywhere  in  the  digit.  Of  particular  interest  is  the 
tip  of  the  digit,  which  undergoes  the  most  extreme  temperature  changes.  Suppose  5°C  is 
selected  as  a  limiting  digit  tip  temperature,  according  to  accepted  experimental  procedures 
involving  cold  exposures  [10].  Under  the  conditions  selected  for  Fig.  3,  it  would  take  about 
44  minutes  for  the  tip  temperature  to  drop  to  5°C.  Another  useful  application  of  the  model 
could  be  the  prediction  of  exposure  times  for  onset  of  frostbite  or  frostnip.  Under  the  present 
set  of  parameters,  Fig.  3  indicates  that  the  minimal  temperature  at  the  digit  tip  would  remain 
above  0°C,  indicating  no  imminent  danger  of  freezing. 

Figures  4  and  5  are  plotted  for  a  case  similar  to  the  previous  one  except  that  both  digit 
base  temperature  and  rate  of  metabolic  heat  generation  are  assumed  to  vary  exponentially 
according  to  Eqs.  (15)  and  (16),  respectively.  One  consequence  of  this  difference  is  in  a 
slightly  larger  time  to  reach  steady  state,  3.28  hrs  compared  to  3.05  hrs.  The  reason  for  this 
prolonged  time  is  due  to  the  lower  final  temperatures  obtained  in  this  case,  Figs.  2  vs.  4. 
Thus,  more  heat  must  be  extracted  and  the  digit  takes  longer  to  reach  a  steady  state. 

Tip  temperatures  of  the  left  little  finger  measured  in  this  laboratory  at  the  nail  bed  under 
sedentary  conditions  [11],  are  compared  to  model  predictions.  The  comparison  is  given  for 
two  different  environmental  temperatures:  -6.7°C  in  Fig.  6  and  0°C  in  Fig.  7.  Of  all 
experimental  conditions  and  subjects  employed  in  this  study  [11],  the  particular  case  selected 
for  comparison  is  the  one  not  involving  any  activity  and  for  the  subject  exhibiting  minimal 
cold-induced  vasodilatation  (CIVD)  response.  This  was  done  since  the  present  model  does 
not  incorporate  this  phenomenon.  It  should  be  pointed  out,  however,  that  the  lack  of  a  CIVD 
response  to  cold  exposure  of  the  fingers  may  be  present  among  a  portion  of  the  population. 
This  makes  such  individuals,  particularly  those  suffering  from  Raynaud's  phenomenon,  more 
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Figure  3.  Temperature  variations  at  different  points  in  the  digit  for  constant  base  temperature 
and  metabolic  heat  generation  rate. 
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and  metabolic  heat  generation  rate. 


vulnerable  to  the  effects  of  cold  stress  [12].  The  lack  of  CIVD  response  is  not,  however, 
limited  to  people  suffering  from  this  disease. 

Parameters  selected  for  calculating  model  predictions  were  based,  as  much  as  was 
possible,  on  actual  values  reported  in  the  study  and  on  established  thermophysical 
properties.  Thus,  a  0.9  clo  (converting  to  h=h,=7.12  W  m  ^C1)  was  used  for  the  thermal 
insulation  applied  on  the  finger  and  finger  tip.  Finger  length  was  set  at  6.5  cm  and  finger 
diameter  at  1 .78  cm.  Initial  finger  base  temperature  was  set  at  33°C  and  was  allowed  to 
drop  to  a  final  value  of  20°C  with  an  exponential  time  constant  of  1.3  hrs,  Eq.  (16). 

The  only  parameter  not  accounted  for  in  the  data  available  was  the  rate  of  metabolic  heat 
generation.  In  the  present  study  this  property  accounts  for  the  combined  effects  of  both  the 
intrinsic  rate  of  metabolic  heat  generation  and  the  amount  of  heat  transported  by  blood 
circulation.  Of  these  two  mechanisms  the  latter  is  by  far  the  dominant  one  and  may  reach 
high  values,  depending  on  the  flow  rate  and  temperature  dif'erence.  Thus,  the  total 
combined  rate  of  heat  generation,  3,  was  estimated  using  minimal  blood  perfusion  rates  of 
0.5  -  1.5  cc/(100  ccmin)  [1,2]  and  a  temperature  difference  of  about  20-25°C.  Based  on 
these  assumptions,  and  the  conformity  of  calculated  to  measured  data,  the  q  used  in  Fig.  6 
was  29,000  W/m3,  decaying  exponentially  to  27,500  W/m3  a  time  constant  of  0.4  hours. 
For  Fig.  7,  the  value  used  was  smaller,  16,500  W which  was  maintained  constant. 

Both  Figs.  6  and  7  show  quite  goo  1  conformity  of  calculated  vs.  measured  data.  The 
largest  deviation,  of  the  order  of  about  1.5°C,  occurs  after  about  30  minutes  in  Fig.  6, 
indicating  fluctuations  in  blood  flow.  In  both  Figs.  6  and  7  there  are  indications  of  weak 
vasomotor  activities  toward  the  end  of  the  experiment.  This  phenomenon  is  particularly 
pronounced  in  Fig.  7  which  otherwise  shows  few  fluctuations  in  finger  tip  temperature. 

The  duration  of  the  two  experiments  was  different,  85  and  105  minutes  in  Figs.  6  and  7, 
respectively.  In  both  experiments,  finger  tip  temperature  dropped  to  about  9°C  at  the  end  of 
the  exposure.  The  model  facilitates  the  estimation  of  the  final  finger  tip  temperatures  and  the 
times  to  reach  these  levels.  Calculations  show  that,  in  both  cases,  the  final  steady  state 
temperatures  would  remain  at  the  safe  level  above  5°C.  Accordingly,  the  finger  tip 
temperatures  would  reach  about  7°C  and  8.4°C  after  about  280  and  237  minutes, 
respectively,  for  the  -6.7°C  and  the  0°C  environments. 

To  further  investigate  endurance  times,  defined  here  as  the  time  for  finger  tip  to  reach 
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Figure  6.  Comparison  of  measured  vs.  calculated  temperature  variations  at  left  little 
finaer  tip  of  subject  #5  at  -6.7°C  environment. 
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Figure  7.  Comparison  of  measured  vs.  calculated  temperature  variations  at  left  little 
finger  tip  of  subject  #y  at  0°C  environment. 


5CC,  a  parametric  analysis  was  performed.  Four  parameters  were  chosen  for  this  analysis: 
neat  transfer  coefficient  of  the  insulation  applied  to  the  digit,  environmental  temperature,  digit 
length  and  digit  diameter.  Results  of  these  analyses  are  shown  in  Figs.  8-10. In  all  cases  a 
5°C  digit  tip  temperature  was  chosen  to  designate  a  limiting  temperature  for  which  an 
endurance  time  may  be  calculated.  In  all  cases  presented  here,  digit  base  temperature  and 
metabolic  heat  generation  rate  varied  exponentially. 

Figures  8-10  indicate  the  trends  for  endurance  times  quite  clearly.  Accordingly,  the  lower 
the  environmental  temperature  (Fig.  8),  the  longer  the  digit  is  (Fig.  9),  the  smaller  the  digit 
diameter  is  (Fig.  10),  the  shorter  the  endurance  times.  An  interesting  result  in  Fig.  9  points 
out  the  fact  that  beyond  a  certain  value  of  the  heat  transfer  coefficient,  endurance  times  are 
no  longer  sensitive  to  digit  length.  In  the  present  case  insulative  values  equal  to  about  0.9 
clo  (h=7.12  W/m2  °C),  or  less,  would  entail  essentially  the  same  endurance  times,  regardless 
of  digit  length. 

In  order  to  evaluate  the  predicted  endurance  times,  two  values  calculated  by  Santee  et  al. 
[11]  are  indicated  by  filled  circles  in  Fig.  8.  These  values  represent  the  average  endurance 
times  obtained  for  the  sedentary  experiments  in  a  -6.7°C  environment  for  two  different 
gloves.  It  is  seen  that  results  of  the  model  underpredict  the  measured  values  by  a  factor  of 
about  2.  This  seemingly  large  disparity  may  easily  be  reconciled  by  recalling  the  fact  that  the 
values  plotted  in  Fig.  8,  were  calculated  for  situations  devoid  of  a  CIVD  response.  The 
measured  values,  however,  were  obtained  for  a  group  of  subjects,  most  of  whom  showed 
active  CIVD  behavior.  Thus,  they  were  able  to  allow  periodic  fluctuations  in  the  blood  supply 
to  the  fingers.  This  increased  the  availability  of  heat  in  the  finger,  considerably  prolonging 
the  endurance  times  over  those  indicated  for  the  non-CIVD  cases. 

Indeed,  when  the  total  heating  rate  is  increased  from  15,000  W/m3,  used  in  Fig.  8,  to 
about  29,000,  used  in  Fig.  6,  model  results  invert  to  overpredict  the  average  measured 
values  by  about  a  factor  of  2.  It  is  to  be  expected,  therefore,  that  the  proper  inclusion  of  the 
CIVD  phenomenon,  would  make  model  predictions  more  realistic  when  applied  to  normal 
populations. 

This  discussion  indicates  the  apparent  limitations  of  the  present  model  and  the  need  for 
further  development.  The  main  usefulness  of  the  present  model  is  that  it  can  provide  a 
systematic  database  for  assessing  lower  bounds  for  endurance  times  under  quite  a  variety  of 
physiological  conditions  as  and  physical  parameters.  These  lower  bounds,  when  calculated 
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Figure  8.  Lower  bounds  of  endurance  times  for  variably  insulated  fingers  exposed  to  different 

environmental  temperatures.  Filled  circles  indicate  average  endurance  times  measured  for 
a  -6.7°C  environment  [11]. 
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Figure  9.  Lower  bounds  of  endurance  times  for  variably  insulated  fingers  of  different  lengths 
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Figure  10.  Lower  bounds  of  endurance  times  for  variably  insulated  fingers  of  different  diameters. 


r 


appropriately,  would  provide  conservative  estimates  of  permissible  exposure  times  of  fingers 
in  cold  environments. 


CONCLUSION 


The  model  developed  in  this  study  is  useful  in  calculating  temperature  variations  along 
insulated  digits  and  in  predicting  their  endurance  times  to  cold  exposures.  Calculated  results 
compared  very  favorably  with  data  obtained  on  human  subjects.  The  present  version  of  the 
model  does  not  include  the  effect  of  cold  induced  vasodilatation.  Expansion  of  the  model  to 
include  this  important  natural  mechanism  will  tend,  as  a  general  rule,  to  predict  longer,  less 
conservative  endurance  times.  Thus,  predictions  based  on  the  present  model  should  be 
useful  whenever  conservative,  lower  bounds  of  endurance  times  to  cold  exposures  are 
required. 
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APPENDIX  A:  DERIVATION  OF  THE  SOLUTION 


For  convenience,  the  set  of  equations  to  be  solved  is  repeated  here: 

±dT  _  *1  -  *H(T-  T)  * 
a  at  dz2  kdK  0  k 

where  a  is  thermal  diffusivity  defined  by  a  =  k/pc. 

T  =  7,(0  @z  =  0,fi0 


(A.1) 


(A.2) 


/*,  (T  -  T0) 


@  z  =  fi,  t  *  0 


(A-3) 


7  =  7,(z)  @  f  =  0,0szsf  (A.4) 

The  following  transformation  is  applied 

r(z,0  =  v(z,0  -  (z)  Fo(0  +  70  (A.5) 

where 

Fo(0  =  7,  (t)  -  T0  (A.6) 

Substitution  of  Eq.  (A.5)  into  Eqs.  (A.1)  -  (A.4)  yields: 
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1  f3v  /  /„x  dFo  ,  _  Sfv  c  Ah, 

a  l¥  “  ^(  )1F  a?  ~  F<,<0  rf?  "  *o‘v  ■  '’(Z)F‘W) 


*(o,o  -  ww  -  /yo  @z  =  o 


*[& -  '  AW'yOl  - 0 


@  z  =  a 


v(z,0)  -  ff(z)  Fo(0)  =  T,(z)  -  T0 


(A.10) 


The  function  f,(z)  may  be  chosen  to  satisfy  certain  conditions  without  loss  of  generality 
[9].  In  the  present  case,  and  for  the  sake  of  convenience,  the  following  conditions  are 
chosen: 


(A.11) 


m  -  -i 


(A.  12) 


k^-  +  MW  =  0 

az 


(A.  13) 


Solution  of  Eqs.  (A.11)  -  (A.  13)  yields 
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k  +  hy  i 


z 


-  1 


(A. 14) 


An  alternative  form  for  f,(z)  would  be 


m 


1  +  Bi  e 


(A.  15) 


where 


Bi  = 


V 

k 


(A.  16) 


and  Bi  is  the  Biot  modulus  indicating  the  ratio  of  convected  heat  at  the  cylinder  tip  to  the  heat 
conducted  toward  the  tip  inside  the  cylinder. 

Substituting  Eq.  (A.  15)  into  Eqs.  (A.7)  -  (A.  10)  and  performing  the  indicated  operations 
yields 


dv 

a~dt 


d?v 

dz2 


tv  *  f,{z)[LFa(t)  . 

a  at 


(A.  17) 


v(0,  t)  =  0  @  z  =  0 


(A.  18) 

| 

i 


k—  *  h.v  =  0  3  z  =  C  (A.19) 

8z 


v(z,0)  =  T,(z)  -  T0  + f,(z)F0(0)  @f  =  0 


(A. 20) 
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Inspection  of  Eqs.  (A.  17)  -  (A.20)  reveals  that  both  boundary  conditions,  Eqs.  (A.  18)  and 
(A.  19),  are  homogeneous.  This  is  the  desired  result  of  the  application  of  the  above 
transformation,  Eq.  (A. 5),  and  the  selection  of  the  conditions  for  f,(z).  One  self-suggesting 
method  of  solution  of  Eqs.  (A.  17)  -(A.20),  would  be  by  applying  an  orthogonal  transformation 
to  eliminate  the  dependence  on  the  spatial  variable,  z.  The  eigenvalue  equation  for  the 
problem  is  obtained  from  Eq.  (A. 19). 

/rp„cosp„  +  VsinPii  =  0  (A.2lj 


or: 


p„cotp„  =  -Bi 


(A.22) 


The  integral  transformation  applied  to  eliminate  the  spatial  variable  z,  is  defined  by  [8,9] 


/ 


ft  ( z ,  t)  K(n,  z)  dz 


(A. 23) 


where  the  bar  indicates  the  transformed  function  which  is  a  function  of  the  remaining  variable 
t  only.  K(n,z)  is  the  kernel  defined  by 


K(n,  z) 


(A. 24) 


where  Kn(z)  are  the  eigenfunctions  of  the  auxiliary  Sturm-Lioville  problem  in  z  [8],  given  for 
the  present  case  by 


C  sin 


P  n* 

l 


(A.25) 
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and  Nn  is  the  normalization  integral 


(A. 26) 


In  the  present  case 


K(n,  z) 


ml  +  Bi2)  9in  m 

+  Bi2  +  Bi)  * 


(A. 27) 


The  transforming  equation  (A.23)  is  applied  now  to  Eqs.  (A.  17)  -  (A.20).  Integration  steps 
would  not  be  repeated  here  but  for  the  sake  of  clarity  the  formal  application  of  the 
transformation  is  shown  in  its  entirety: 


{  dz 2 


~  f  v(z,  t)  K{n,  z)dz  =  {  ^  K{n,  z)dz  -  if  v  K(n,  z)dz 


dF  * 

/  K[n.z)dz  ♦  {LFa{t)  .  ±-^)  f/,MK(n.z)dz 


(A.28) 


Both  boundary  conditions,  Eqs.  (A. 18)  and  (A.  19),  are  satisfied  by  the  set  of 
eigenfunctions.  The  initial  condition,  Eq.  (A.20)  is  transformed  by  applying  Eq.  (A.23)  to  yield 

fv(z,0)K[n,z)dz  =  \(T,{z)  -  T0  +  f^z)Fo(0))K(n,z)dz  (A.29) 

0  0 

Performing  the  indicated  integrations  and  rearranging  terms,  yields  an  ordinary  differential 
equation  of  the  transformed  function  v(t): 
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-aG(t) 


(A. 30) 


dv 

dt 


ak2n\ 


subject  to  the  transformed  initial  condition 


v(0)  =  T  @  t  =  0 


(A-31) 


where 


o 


(A.32) 


<3(0  -  . 


1  dFol  <*f)C0„ 


dt 


*P, 


(cosp„  -  1 ) 


(A.33) 


2  (Pi  +  Bi2) 

«“(P 5  -  -  a) 


(A.34) 


and 


T  =  S„/ [  r/z)  -  rjsin^Jz  -  i^FJO) 


(A.35) 


The  solution  of  Eq.  (A.30),  subject  to  Eq.  (A. 31),  is  straightforward  and  is  given  by 
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(A. 36) 


t 

v(t)  «  (x  -  a  j  G(ri)exp(a  A.zri)dri  }exp(-o>,2f) 

rj  *=0 


Inversion  of  the  transformed  function,  v(t),  is  achieved  by  employing  the  orthogonality 
property  of  the  eigenfunctions.  This  inversion  is  performed  by 


v  (z,  t)  -  £  v(t)Kn(z) 

n*  1 


(A.37) 


v(*.  0 


fEh 


D 

«/  C?(r) )  exp  (a  A,2ti)  dii)  exp(-a  A.2f)sin-y 


(A.38) 


The  solution  for  the  temperature  T(z,t)  is  obtained  by  substituting  Eq.  (A.38)  into  Eq. 
(A.5),  to  yield 


T(z,  t)  =  (Eli  -  a  fG(ii)exp(a  J.2Ti)dn  )exp(-A.2a0sin 

n-1  * 


M 


(A.39) 


-  Ta  -  m\TM  -  ro] 


Equation  (A.39)  is  the  solution  for  the  original  set  of  equations  (A.1)  -  (A. 4). 


35 


APPENDIX  B:  COMPUTER  PROGRAM  PRINTOUT  AND  SAMPLE  OUTPUT 


PROGRAM  TEMP4;  {  Calculation  of  the  time-dependent  temperature  field 
in  a  cylinder  in  the  axial  direction  for  variable  (exponential) 
temperature  at  z=0  and  variable  (exponential)  heat  generation. 
Program  written  in  Turbo-Pascal  Version  5.5  } 

{  Current  version  employs  time  in  hours  } 

uses  crt, printer; 

LABEL 

30,40,50, 

100,200,300,400,600; 

CONST  {  definition  of  parameters  } 

p-3.1 4159265; 

k=0.418;  { thermal  conductivity,  W/m/C  } 

h=2.0;  {  heat  transfer  coefficient  on  circumference,  W/mA2/C  } 

hi  =2.0;  {  heat  transfer  coefficient  on  tip,  W/mA2/C  } 

D=0.02;  {  cylinder  diameter,  m  } 

qi=2000.;  { initial  heat  generation  (metabolic)  rate,  W/mA3  } 
qf=200;  { final  heat  generation  (metabolic)  rate,  W/mA3  } 

T0=-5.0;  {  environmental  temperature,  C  } 

T 1  =30. ;  {  constant  base  temperature,  C  } 

T2=30.0;  {  base  temperature,  linear  distribution,  C  } 

T3=20.0;  { tip  temperature,  linear  distribution,  C  } 

1=0.12;  {  cylinder  length,  m  } 

M=40;  {  number  of  terms  in  series  (maximum) } 

ALF=0. 0004546;  { thermal  diffusivity,  mA2/hr } 
delta=1 .3;  { time  constant  for  base  temperature,  hr } 

Tfin=20.0;  { final  base  temperature,  C  } 

Tinit=30.0;{  initial  base  temperature,  C  } 

eps=1 .3;  { time  constant  for  metabolic  heat  rate,  hr } 

KK-  11;  {  number  of  time  steps  in  minutes  } 

TYPE 

vec=ARRAY[1  ..40]  OF  double; 
mat=ARRAY[1  ..20,1  ..41]  OF  double; 


VAR 

ft,ct  .text; 

Tfile.cfile  :string; 

init, uniform, linear  :string; 

l,J,N  integer; 

T  :mat; 

Bn.Gn.Hn  :vec; 
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Cn,LAMSn,TAUn,COEFn,b  :vec; 

TIME, temp, Tmax  :double; 

Bi,con,asum,ax  :double; 

x,c,LI,y  :double; 

del,ddx,dev,u,uu,dx  idouble; 

q,qq,f,g,A,tot  :doubie; 

2,ff1,SUM,SUM1,SU,fsin  idouble; 

FUNCTION  Par(y,c:double):double; 
begin 

Par:=y*COS(y)/SIN(y)+c 

end; 

FUNCTION  Der(y:double):double; 
begin 

Der:=COS(y)/SIN(y)-y*(1  .-sqr(COS(y)/SIN(y))) 
end; 

FUNCTION  F1(y:double):double; 
begin 

FI  :=Bi/(1.+Bi)*(y/l)-1. 
end; 

BEGIN  {  begin  execution  of  main  program  } 

ClrScr; 

init:=’linear’; 

{  Nullify  initial  values  of  vectors  } 

FOR  n:=1  TO  40  DO 
BEGIN 
Bn[n]  :=0.; 

Cn[n]:=0.; 

LAMSn[n]:=0.; 

TAUn[n]:=0.; 

COEFn[n]:=0.; 

b[n]:=0.; 

Gn[n]:=0.; 

Hn[n]:=0. 

END; 

WRITELN('k=  ’,k:6:2,’  h=*,h:6:2,  *  h1=’,h1 :6:2,'  ALF=’,ALF:16:8) 

WRITELN(’qi=  \qi:6:2,’  qf=’,qf:6,’  T0=’,T0:6:2); 

WRITELN(’T1=  ’,T1:6:2,’  T2=  ’,T2:6:2,'  T3=’,T3:6:2)  ; 
WRITELN(’I=  ’,1:6:2,’  D=’,D:6:2)  ; 


{  Procedure  to  calculate  the  eigenvalues  of: 
f(x)=x*cot(x)+c } 

Bi:=hri/k;  {  Biot  modulus  } 

LI:=4.*h/k/D; 

c:=Bi; 

WRITELN(’Biot=’,c:6:2)  ; 
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WRITELN  ; 

{  Check  values  of  the  parameter  c  and  set  initial  guesses  and  increments  } 
IF  c<1 .37 

THEN  GOTO  50 

ELSE 

begin 

x:=2.2; 

ddx:=2.; 

del:=0.0001*c 

end; 

IF  0=1.37  THEN 
begin 
x:=2.1 ; 
ddx:=2. 
end  ; 

IF  c>=50.  THEN 
begin 
x:=2.7; 
ddx:=3. 
end; 

IF  c>=200.  THEN 
begin 
x:=3.; 
ddx:=3.1 
end; 

IF  c>=2000.  THEN 
begin 
x:=3.  ; 
ddx:=3.14 
end  ; 


{  calculate  the  first  iterative  values  } 

FOR  i:=1  TO  M  DO 
BEGIN 
dev:=1.  ; 

30:  u:=1 .  ; 
uu:=1.  ; 
dx:=0. 01/dev; 

WHILE  (u*uu)>0  DO 
begin 
x:=x+dx; 
u:=Par(x,c)  ; 
uu:=Par(x+dx,c) 
end  ; 

IF  ABS(u)<del  THEN 
ELSE 
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begin 

dev:=dev*100.  ; 

IF  dev>1E08  THEN  GOTO  600; 

GOTO  30 
end  ; 

600:  b[i]:=x  ; 
x:=x+ddx 
END  ; 

GOTO  40  ; 

{  Calculate  approximate  values  of  roots  by  using  the 
slope  of  cot(x)  at  (2*n-1)*p/2  } 

50:  FOR  n:=1  TO  M  DO  {  calculation  for  low  values  of  c  } 

BEGIN 

IF  (n=1)  AND  (c<(0.5-1 ))  THEN 

b[1]:=0.1 

ELSE 

{  This  portion  calculates  the  eigenvalues  for  low  and  negative  values 
of  the  variable  } 

begin 

q:=(2*n-1)*pi/2  ; 
qq:=q*q+4.*c  ; 

b[n]:=(q+SQRT(qq))/2. 
end 
END; 

{  Start  Newton’s  procedure  } 

FOR  i:=1  TO  M  DO 

BEGIN 

f:=1  • ; 

x:=b[i]  ; 

del:=0.0001  *EXP(INT(0.1*i)*LN{1 0))  ; 

WHILE  ABS(f)>del  DO 
begin 

f:=Par(x,c)  ; 
g:=Der(x)  ; 
x:=x-f/g 
end  ; 
b[i]:=x+f/g 
END; 

40:  {FOR  l:=1  TO  M  DO 
WRITEIn{l,b[l],Par(b[l],c));} 

{  End  of  procedure  for  calculating  eigenvalues  } 

FOR  l:=1  TO  M  DO 
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BEGIN  {  calculation  of  parameters  in  series  solution  } 

A:=b[l]*b[l]+Bi*Bi; 

Bn[l]:=2*A/(IT(A+Bi))  ; 

Cn[l]:=Bn[l]*l/b[l]*(qf*(COS(b[l])-1 .0)/k+U*(Tfin-T0))  ; 
LAMSn[l]:=b[l]*b[l]/l/l+LI  ; 

TAUn[l]:=Bn[l]*l/b[l]*((T3-T2)*SIN(b[l])/b[l]+(T2-Tinit)-(T3-TO)*COS(b[l]»; 

COEFn[l]:=Cn[l]+LAMSn[l]*TAUn[l]; 

Gn[l]:=Bn[l]*l/b[l]*LAMSn[l]*(Tinit-Tfin)*(ALF*DELTA*LI-1.0)/(ALF*LAMSn[l]* 

DELTA-1.0); 

Hn[l]:=ALF*l*Bn[l]/k/b[l]'(COS(b[l])-1.0)'LAMSn[l]*(qi-qf)*eps/ 
(eps*ALF*LAMSn[l]-1 .0) 

END; 

{  estimation  of  time  to  reach  steady  state  } 

Tmax:=LN(abs(1000*(COEFn[1]+Gn[1]+Hn[1])/Cn[1]))/(ALF*LAMSn[1]); 
WRITELNf  Tmax=’1Tmax:10:5,’  hr',’  T0=’,T0:10:5,’  deg  C’); 

Tfile  :=  ’Test.dat’;  {  open  output  file  } 

assign(ft.Tfile); 

rewrite(ft); 

WRITELN(ft,'k=  ’,k:6;2,’  h=’,h:6:2,  ’  h1=',h1 :6:2,’  ALF=’,ALF:16:8)  ; 

WRITELN(ft,’qi=  ',qi:6:2,'  qf=’,qf:6>’  T0=‘,T0:6:2); 

WRITELN(ft,’T1=  ’,T1;6;2,’  T2=  *,T2:6;2>’  T3=’,T3:6;2)  ; 

WRITELN(ft,’l= ’,1:6:2,’  D=’,D:6:2)  ; 

WRITELN(ft,’Biot=’,c:6:2)  ; 

WRITELN(ft); 

FOR  i:=1  TO  KK  DO  { time  steps  } 
begin 

write(ft,i-1:2,’  ’); 

TIME:=(i-1); 

IF  TIME=0.  THEN 
TIME:=0. 00001 
ELSE; 

if  i=KK  then  TIME>Tmax 
else  ; 

FOR  j:=1  TO  1 1  DO 
BEGIN 
z:=(j-1  )*l/10.  ; 
ff  1  :=F1  (z)  ; 

SUM:=0.  ; 

con:=1.  ; 

N:=1  ; 

WHILE  ((N<(M+1))  AND  (con>asum))  DO 
begin  {while} 

SUM1:=SUM  ; 

IF  (abs(ALF*LAMSn[N]‘TIME)>LN{abs(1000*(COEFn[N]+Gn[N]+Hn[n]) 

/Cn[N])))  THEN 
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SUM:=SUM-(Cn[N]+Gn[N]*EXP(-TIME/DELTA)+Hn[N]*EXP(-TIME/eps))* 

SIN(b[N]*z/l)/LAMSn[N] 

ELSE 

SUM:=SUM+((COEFn[N]+Gn[N]+Hn[N])*EXP(-ALF‘LAMSn[N]*TIME)-Cn[N]-Gn[N]* 
EXP(-TIME/DELTA)-Hn[N]*EXP(-TIME/eps))*SIN(b[N]*z/l)/LAMSn[N]; 
con:=abs(ABS(SUM)-ABS(SUM1 ))  ; 
asum:=0.001  *ABS(SUM)  ; 

N:=N+1 


end  ;  {while} 

200:  {WRITE(z,N,SUM);} 

TEMP:=l‘SUM+T0-ff1  *(Tfin+(Tinit-Tfin)*EXP(-TIME/DELTA)-TO)  ; 

T[j,i]:=TEMP; 

write(ft,  temp:6:4,’  ’); 

END  ; 
writeln(ft); 

end;  {  end  time  steps  } 

close(ft); 


tot:=0.; 

for  i:=1  to  KK  do  { time  steps  } 
begin 

write(i-1 :2); 

for  j:=1  to  1 1  do  {  space  steps  } 
begin 

write  (T[j ,  i]  :7 :3) ; 

tot:=tot+T[j,i]  {  sum  of  all  printed  values;  used  to  check  consistency  } 

end; 

writeln 

end;  {  end  printing  of  temperature  matrix  } 
writeln(’  totals’, tot:1 0:5)  ; 


100:  END. 
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Sample  output  of  PROGRAM  TEMP4 
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20.0002  12.3730  7.2106  3.6390  1.0699  -0.6014  -1.7664  -2.5289  -3.0415  -3.3643  -3.5076 


GLOSSARY 


A  -  cross  sectional  area 

Bi  -  Biot  modulus,  Eq.  (A.  16) 

Bn  -  parameter,  Eq.  (7) 

c  -  specific  heat 

Cn  -  parameter,  Eq.  (18) 

d  -  diameter 

f,  -  function,  Eq.  (10) 

F0  -  function,  Eq.  (9) 

G  -  defined  by  Eq.  (8) 

Gn  -  parameter,  Eq.  (19) 

h  -  heat  transfer  coefficient  at  cylinder  circumference 

h,  -  heat  transfer  coefficient  at  cylinder  tip 

Hn  -  eigenfunction,  Eq.  (A.25) 

k  -  thermal  conductivity 

K  -  kernel,  Eg.  (A.24) 

Kn  -  eigenfunction,  Eq.  (A.25) 

I  -  cylinder  length 

L  -  4h/kd 

n  -  integer 

Nn  -  normalization  integral,  Eq.  (A.26) 

p  -  cylinder  circumference 
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q 


-  metabolic  heat  generation  rate  per  unit  volume 

-  time 


t 

T  -  temperature 

T2,  T3  -  temperatures  in  linear  initial  distribution,  Eq.  (13) 

v  -  function,  Eq.  (A. 5) 

z  -  axial  coordinate 


GREEK  LETTERS 

a  -  k/pc,  thermal  diffusivity 

Pn  -  eigenvalue,  Eq.  (1 1) 

5  -  time  constant  for  cylinder  base  temperature,  Eq.  (15) 

e  -  time  constant  for  metabolic  heat  generation  rate,  Eq.  (16) 

C,  -  integration  variable,  Eq.  (12) 

T|  -  integration  variable,  Eq.  (5) 

-  parameter,  Eq.  (A.32) 

p  -  density 

x  -  defined  by  Eq.  (6) 

Q  -  fjnction 

Qn  -  transformed  function,  Eq.  (A. 23) 

SUBSCRIPTS 


o  -  environmental 

1  -  cylinder  base 
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f  -  final 

i  -  initial 


SUPERSCRIPTS 


-  transformed  function,  defined  by  Eq.  (A.23) 

-  steady  state 
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